# Setup -------------------------------------------------------------------
rm(list = ls())
source("figure_design.R")

# Load data ---------------------------------------------------------------

data <- readRDS("WVS_turkey_v1_6_subset.rds") %>%
  mutate_at(vars(Q224:Q234), ~5 - .) %>%
  mutate(J_INTDATE = as.Date(as.character(J_INTDATE), format = "%Y%m%d"),
         running = difftime(J_INTDATE, time2 = as.Date("2018-04-18"), units = "days"),
         running = as.numeric(running, units = "days"),
         D_INTERVIEW = running) %>%
  filter(D_INTERVIEW != 0) %>%
  mutate(D_INTERVIEW = ifelse(D_INTERVIEW > 0, D_INTERVIEW - 1, D_INTERVIEW),
         treated = ifelse(D_INTERVIEW >= 0, 1, 0)) %>%
  mutate(H_URBRURAL = car::recode(H_URBRURAL, "1=0;2=1"),
         Q260 = as.factor(Q260),               # Gender
         Q279 = car::recode(Q279, "7=1;else=0"),  # Employment status
         Q202 = car::recode(Q202, "1=1;else=0"))  # TV

# Subset ------------------------------------------------------------------

data %>% filter(Partyabb == "AKP"|Partyabb == "MHP") -> data_incumbent
data %>% filter(Partyabb== "SHP/CHP"|Partyabb == "HDP") -> data_opposition

# Figure 5a ---------------------------------------------------------------

create_lm <- function(data, data_incumbent, data_opposition, dep_var) {
  lm1 <- lm(as.formula(paste0(dep_var, " ~ treated * D_INTERVIEW")), data = data[data$D_INTERVIEW >= -4 & data$D_INTERVIEW <= 4, ])
  lm2 <- lm(as.formula(paste0(dep_var, " ~ treated * D_INTERVIEW")), data = data_incumbent[data_incumbent$D_INTERVIEW >= -4 & data_incumbent$D_INTERVIEW <= 4, ])
  lm3 <- lm(as.formula(paste0(dep_var, " ~ treated * D_INTERVIEW")), data = data_opposition[data_opposition$D_INTERVIEW >= -5 & data_opposition$D_INTERVIEW <= 5, ])
  return(list(lm1, lm2, lm3))
}

# List of dependent variables
dependent_vars <- paste0("Q", 224:234)

# List to store the linear models
lm_list <- lapply(dependent_vars, function(dep_var) create_lm(data, data_incumbent, data_opposition, dep_var))

# Define the confidence levels
confidence_levels <- c(0.90, 0.95)

# Initialize the dataframe
plotdata <- data.frame(
  y = rep(1:11, each = 2),
  Group = rep(c("Incumbent supporters", "Opposition supporters"), times = 11),
  x = rep(1:11, each = 2),
  xmin95 = NA,
  xmax95 = NA,
  xmin90 = NA,
  xmax90 = NA,
  question = rep(paste0("Q", 224:234), each = 2),
  label = rep(c(
    "Votes counted fairly",
    "Opposition prevented from running",
    "TV favors governing party",
    "Voters bribed",
    "Journalists provide fair coverage",
    "Election officials fair",
    "Rich people buy elections",
    "Voters threatened with violence",
    "Voters offered genuine choice",
    "Women have equal opportunities to run",
    "Honest elections make a difference"
  ), each = 2)
)

for (i in seq_along(dependent_vars)) {
  dep_var <- dependent_vars[i]
  lm_model_incumbent <- lm_list[[i]][[2]]  # Access the second linear model
  lm_model_opposition <- lm_list[[i]][[3]]  # Access the third linear model
  
  # Store the coefficients for incumbent supporters
  plotdata$x[i * 2 - 1] <- coef(lm_model_incumbent)["treated"]
  
  # Store the coefficients for opposition supporters
  plotdata$x[i * 2] <- coef(lm_model_opposition)["treated"]
  
  # Store the confidence intervals for incumbent supporters
  conf_int_90_incumbent <- confint(lm_model_incumbent, level = confidence_levels[1])
  plotdata$xmin90[i * 2 - 1] <- conf_int_90_incumbent["treated", 1]
  plotdata$xmax90[i * 2 - 1] <- conf_int_90_incumbent["treated", 2]
  
  conf_int_95_incumbent <- confint(lm_model_incumbent, level = confidence_levels[2])
  plotdata$xmin95[i * 2 - 1] <- conf_int_95_incumbent["treated", 1]
  plotdata$xmax95[i * 2 - 1] <- conf_int_95_incumbent["treated", 2]
  
  # Store the confidence intervals for opposition supporters
  conf_int_90_opposition <- confint(lm_model_opposition, level = confidence_levels[1])
  plotdata$xmin90[i * 2] <- conf_int_90_opposition["treated", 1]
  plotdata$xmax90[i * 2] <- conf_int_90_opposition["treated", 2]
  
  conf_int_95_opposition <- confint(lm_model_opposition, level = confidence_levels[2])
  plotdata$xmin95[i * 2] <- conf_int_95_opposition["treated", 1]
  plotdata$xmax95[i * 2] <- conf_int_95_opposition["treated", 2]
}

coeffplot <- ggplot(plotdata,aes(x = label, y = x, color = Group, shape = Group)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_linerange(aes(ymin = xmin90,
                     ymax = xmax90),
                 linewidth = 1,
                 position = position_dodge(0.4)) +
  geom_linerange(aes(ymin = xmin95,
                     ymax = xmax95),
                 linewidth = 0.45,
                 position = position_dodge(0.4)) + 
  geom_point(size=3,position = position_dodge(0.4)) + 
  xlab("") +
  ylab("Effect of Announcement (LATE)") +
  theme(legend.title=element_blank()) +
  coord_flip() + 
  theme(axis.text.x = element_text(size = 8, angle = 45,hjust = 1, vjust = 1)) + 
  scale_color_manual(values = plot_colors)  +
  theme(legend.position = "top",
        text = element_text(size = 15),
        axis.text = element_text(size = 14))

pdf(file = "fig_5a.pdf",width = 7.5,height = 6,onefile=F)

coeffplot

dev.off()

# Figure 5b ---------------------------------------------------------------

# Function to run rdrobust and retrieve results
run_rd_robust <- function(data, data_incumbent, data_opposition, running, running_incumbent, running_opposition) {
  rdrb_full_95 <- rdrobust(data, running, c = 0, masspoints = "adjust", all = TRUE, level = 95)
  rdrb_full_90 <- rdrobust(data, running, c = 0, masspoints = "adjust", all = TRUE, level = 90)
  rdrb_in_95 <- rdrobust(data_incumbent, running_incumbent, c = 0, masspoints = "adjust", all = TRUE, level = 95)
  rdrb_in_90 <- rdrobust(data_incumbent, running_incumbent, c = 0, masspoints = "adjust", all = TRUE, level = 90)
  rdrb_op_95 <- rdrobust(data_opposition, running_opposition, c = 0, masspoints = "adjust", all = TRUE, level = 95)
  rdrb_op_90 <- rdrobust(data_opposition, running_opposition, c = 0, masspoints = "adjust", all = TRUE, level = 90)
  result <- list(rdrb_full_95, rdrb_full_90, rdrb_in_95, rdrb_in_90, rdrb_op_95, rdrb_op_90)
  return(result)
}

# List of dependent variables
dependent_vars <- paste0("Q", 224:234)

# Run rdrobust for each dependent variable
rd_results <- lapply(
  dependent_vars,
  function(dep_var) {
    run_rd_robust(
      data[[dep_var]], data_incumbent[[dep_var]], data_opposition[[dep_var]],
      data$D_INTERVIEW, data_incumbent$D_INTERVIEW, data_opposition$D_INTERVIEW
    )
  }
)

# Initialize the dataframe
plotdata_robust <- data.frame(
  y = rep(1:11, each = 2),
  Group = rep(c("Incumbent supporters", "Opposition supporters"), times = 11),
  x = rep(1:11, each = 2),
  xmin95 = NA,
  xmax95 = NA,
  xmin90 = NA,
  xmax90 = NA,
  question = rep(paste0("Q", 224:234), each = 2),
  label = rep(c(
    "Votes counted fairly",
    "Opposition prevented from running",
    "TV favors governing party",
    "Voters bribed",
    "Journalists provide fair coverage",
    "Election officials fair",
    "Rich people buy elections",
    "Voters threatened with violence",
    "Voters offered genuine choice",
    "Women have equal opportunities to run",
    "Honest elections make a difference"
  ), each = 2)
)

# Populate the dataframe with coefficients and confidence intervals
for (i in seq_along(dependent_vars)) {
  plotdata_robust$x[i * 2 - 1] <- rd_results[[i]][[3]]$coef[3]
  plotdata_robust$x[i * 2] <- rd_results[[i]][[5]]$coef[3]
  
  plotdata_robust$xmin90[i * 2 - 1] <- rd_results[[i]][[4]]$ci["Robust", "CI Lower"]
  plotdata_robust$xmax90[i * 2 - 1] <- rd_results[[i]][[4]]$ci["Robust", "CI Upper"]
  plotdata_robust$xmin95[i * 2 - 1] <- rd_results[[i]][[3]]$ci["Robust", "CI Lower"]
  plotdata_robust$xmax95[i * 2 - 1] <- rd_results[[i]][[3]]$ci["Robust", "CI Upper"]
  
  plotdata_robust$xmin90[i * 2] <- rd_results[[i]][[6]]$ci["Robust", "CI Lower"]
  plotdata_robust$xmax90[i * 2] <- rd_results[[i]][[6]]$ci["Robust", "CI Upper"]
  plotdata_robust$xmin95[i * 2] <- rd_results[[i]][[5]]$ci["Robust", "CI Lower"]
  plotdata_robust$xmax95[i * 2] <- rd_results[[i]][[5]]$ci["Robust", "CI Upper"]
}

coeffplot <- ggplot(plotdata_robust,aes(x = label, y = x, color = Group, shape = Group)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_linerange(aes(ymin = xmin90,
                     ymax = xmax90),
                 linewidth = 1,
                 position = position_dodge(0.4)) +
  geom_linerange(aes(ymin = xmin95,
                     ymax = xmax95),
                 linewidth = 0.45,
                 position = position_dodge(0.4)) + 
  geom_point(size=3,position = position_dodge(0.4)) + 
  xlab("") +
  ylab("Effect of Announcement (LATE)") +
  theme(legend.title=element_blank()) +
  coord_flip() + 
  theme(axis.text.x = element_text(size = 8, angle = 45,hjust = 1, vjust = 1)) + 
  scale_color_manual(values = plot_colors) +
  theme(legend.position = "top",
        text = element_text(size = 15),
        axis.text = element_text(size = 14))

  coeffplot

pdf(file = "fig_5b.pdf",width = 7.5,height = 6,onefile=F)

coeffplot

dev.off()

# Appendix C --------------------------------------------------------------

# Table C.2 ---------------------------------------------------------------

data %>% 
  select(Q224,Q225,Q226,Q227,Q228,Q229,Q230,Q231,Q232,Q233,Q234,
         Q262,Q260,H_URBRURAL,Q279,Q202) %>%
  mutate_if(is.factor, as.numeric) -> 
  data_descriptives

stargazer(data_descriptives, covariate.labels = c("Votes counted fairly",
                                                  "Opposition prevented from running",
                                                  "TV favors governing party",
                                                  "Voters bribed",
                                                  "Journalists provide fair coverage",
                                                  "Election officials fair",
                                                  "Rich people buy elections",
                                                  "Voters threatened with violence",
                                                  "Voters offered genuine choice",
                                                  "Women have equal opportunities to run",
                                                  "Honest elections make a difference",
                                                  "Age",
                                                  "Gender",
                                                  "Urban/rural",
                                                  "Unemployment",
                                                  "TV Consumption"),
          out = "tab_c2.tex")

# Figure C.1 --------------------------------------------------------------

# List of plot titles and corresponding column names
plot_info <- list(
  list(title = "Votes are counted fairly", col_name = "Q224"),
  list(title = "Opposition candidates are prevented from running", col_name = "Q225"),
  list(title = "TV news favors the governing party", col_name = "Q226"),
  list(title = "Voters are bribed", col_name = "Q227"),
  list(title = "Journalists provide fair coverage of elections", col_name = "Q228"),
  list(title = "Election officials are fair", col_name = "Q229"),
  list(title = "Rich people buy elections", col_name = "Q230"),
  list(title = "Voters are threatened with violence at the polls", col_name = "Q231"),
  list(title = "Voters are offered a genuine choice in the elections", col_name = "Q232"),
  list(title = "Women have equal opportunities to run the office", col_name = "Q233"),
  list(title = "Honest elections make a difference",col_name = "Q234")
)

# Loop over the list of plot information
for (i in seq_along(plot_info)) {
  # Open a PDF device for each plot
  pdf(paste0("fig_c1", letters[i], ".pdf"), width = 8, height = 4)
  
  # Extract data for the current plot
  current_data <- data_incumbent[data_incumbent$D_INTERVIEW >= -10 & data_incumbent$D_INTERVIEW <= 10, ]
  
  # Plot the data
  rdplot(current_data[[plot_info[[i]]$col_name]], current_data$D_INTERVIEW, ci = 95, c = 0,
         title = "", y.label = plot_info[[i]]$title, x.label = "Days +/- election announcement")
  
  # Close the current plot device
  dev.off()
}

# Figure C.2 --------------------------------------------------------------

# List of plot titles and corresponding column names
plot_info <- list(
  list(title = "Votes are counted fairly", col_name = "Q224"),
  list(title = "Opposition candidates are prevented from running", col_name = "Q225"),
  list(title = "TV news favors the governing party", col_name = "Q226"),
  list(title = "Voters are bribed", col_name = "Q227"),
  list(title = "Journalists provide fair coverage of elections", col_name = "Q228"),
  list(title = "Election officials are fair", col_name = "Q229"),
  list(title = "Rich people buy elections", col_name = "Q230"),
  list(title = "Voters are threatened with violence at the polls", col_name = "Q231"),
  list(title = "Voters are offered a genuine choice in the elections", col_name = "Q232"),
  list(title = "Women have equal opportunities to run the office", col_name = "Q233"),
  list(title = "Honest elections make a difference",col_name = "Q234")
)


# Loop over the list of plot information
for (i in seq_along(plot_info)) {
  # Open a PDF device for each plot
  pdf(paste0("fig_c2", letters[i], ".pdf"), width = 8, height = 4)
  
  # Extract data for the current plot
  current_data <- data_incumbent[data_opposition$D_INTERVIEW >= -10 & data_incumbent$D_INTERVIEW <= 10, ]
  
  # Plot the data
  rdplot(current_data[[plot_info[[i]]$col_name]], current_data$D_INTERVIEW, ci = 95, c = 0,
         title = "", y.label = plot_info[[i]]$title, x.label = "Days +/- election announcement")
  
  # Close the current plot device
  dev.off()
}


# Tables C.3 - C.13 -------------------------------------------------------

# Define output filenames
output_files <- character(13)
for (i in 3:13) {
  output_files[i-2] <- paste0("tab_c", i, ".tex")
}

# Loop over lm_list and generate corresponding tables
for (i in seq_along(lm_list)) {
  texreg(lm_list[[i]],
            file = output_files[i], # Specify output filename
            include.ci = FALSE,
            digits = 3,
            custom.coef.map = list(
              "(Intercept)" = "Intercept",
              "treated" = "Treated",
              "D_INTERVIEW" = "Running variable",
              "treated:D_INTERVIEW" = "Treated $\\times$ running variable"
            ),
            custom.header = list("Pooled" = 1,
                                 "Incumbent" = 2,
                                 "Opposition" = 3),
            booktabs = TRUE,
            use.packages = FALSE,
            include.rmse = FALSE
  )
}

# Table C.14 --------------------------------------------------------------

items <- c(
  "Votes counted fairly",
  "Opposition prevented from running",
  "TV favors governing party",
  "Voters bribed",
  "Journalists provide fair coverage",
  "Election officials fair",
  "Rich people buy elections",
  "Voters threatened with violence",
  "Voters offered genuine choice",
  "Women have equal opportunities to run",
  "Honest elections make a difference"
)

# Extract estimates, standard errors, p-values, and bandwidths
estimates <- sapply(rd_results[1:11], function(x) x[[3]][["coef"]]["Robust", "Coeff"])
standard_errors <- sapply(rd_results[c(1:11)], function(x) x[[3]][["se"]]["Robust", "Std. Err."])
p_values <- sapply(rd_results[1:11], function(x) x[[3]][["pv"]]["Robust", "P>|z|"])
bandwidths <- sapply(rd_results[1:11], function(x) x[[3]][["bws"]]["h","left"])

# Create data frame
df_robust_estimates <- data.frame(
  item = rep(items),
  estimate = unlist(estimates),
  SE = unlist(standard_errors),
  pvalue = unlist(p_values),
  bandwidth = unlist(bandwidths)
)

write.csv(df_robust_estimates,file = "tab_c14.csv",row.names = F)

# Table C.15 --------------------------------------------------------------

items <- c(
  "Votes counted fairly",
  "Opposition prevented from running",
  "TV favors governing party",
  "Voters bribed",
  "Journalists provide fair coverage",
  "Election officials fair",
  "Rich people buy elections",
  "Voters threatened with violence",
  "Voters offered genuine choice",
  "Women have equal opportunities to run",
  "Honest elections make a difference"
)

# Extract estimates, standard errors, p-values, and bandwidths
estimates <- sapply(rd_results[1:11], function(x) x[[5]][["coef"]]["Robust", "Coeff"])
standard_errors <- sapply(rd_results[c(1:11)], function(x) x[[5]][["se"]]["Robust", "Std. Err."])
p_values <- sapply(rd_results[1:11], function(x) x[[5]][["pv"]]["Robust", "P>|z|"])
bandwidths <- sapply(rd_results[1:11], function(x) x[[3]][["bws"]]["h","left"])

# Create data frame
df_robust_estimates <- data.frame(
  item = rep(items),
  estimate = unlist(estimates),
  SE = unlist(standard_errors),
  pvalue = unlist(p_values),
  bandwidth = unlist(bandwidths)
)

write.csv(df_robust_estimates,file = "tab_c15.csv",row.names = F)

# Figure C.3 --------------------------------------------------------------

rdd <- rdd_data(y = data_incumbent$Q232, x= data_incumbent$D_INTERVIEW, cutpoint=0)
reg_np <- rdd_reg_np(rdd, inference = "np",bw = 5) #Dataprep

pdf(file = "fig_c3.pdf",
    width = 6,
    height = 5,
    onefile = F)

plotSensi(reg_np, from = 3, to = 15) # BW Sensitivity

dev.off()

# Figure C.4a -------------------------------------------------------------

create_lm <- function(data, data_incumbent, data_opposition, dep_var) {
  lm1 <- lm(as.formula(paste0(dep_var, " ~ treated * D_INTERVIEW")), data = data[data$D_INTERVIEW >= -4 & data$D_INTERVIEW <= 4, ])
  lm2 <- lm(as.formula(paste0(dep_var, " ~ treated * D_INTERVIEW")), data = data_incumbent[data_incumbent$D_INTERVIEW >= -4 & data_incumbent$D_INTERVIEW <= 4, ])
  lm3 <- lm(as.formula(paste0(dep_var, " ~ treated * D_INTERVIEW")), data = data_opposition[data_opposition$D_INTERVIEW >= -5 & data_opposition$D_INTERVIEW <= 5, ])
  return(list(lm1, lm2, lm3))
}

# List of dependent variables
dependent_vars <- c("Q49","Q59","Q142","Q147","Q170","Q199","Q222")

# List to store the linear models
lm_list <- lapply(dependent_vars, function(dep_var) create_lm(data, data_incumbent, data_opposition, dep_var))

# Define the confidence levels
confidence_levels <- c(0.90, 0.95)

# Initialize the dataframe
plotdata <- data.frame(
  y = rep(1:7, each = 2),
  Group = rep(c("Incumbent supporters", "Opposition supporters")),
  x = rep(1:7, each = 2),
  xmin95 = NA,
  xmax95 = NA,
  xmin90 = NA,
  xmax90 = NA,
  question = rep(c(dependent_vars), each = 2),
  label = rep(c(
    "Satisfaction with life",
    "Trust in neighborhood",
    "Fear of losing job",
    "Worried about terrorism",
    "Only own religion acceptable",
    "Interest in politics",
    "Frequency of voting"
  ), each = 2)
)

for (i in seq_along(dependent_vars)) {
  dep_var <- dependent_vars[i]
  lm_model_incumbent <- lm_list[[i]][[2]]  # Access the second linear model
  lm_model_opposition <- lm_list[[i]][[3]]  # Access the third linear model
  
  # Store the coefficients for incumbent supporters
  plotdata$x[i * 2 - 1] <- coef(lm_model_incumbent)["treated"]
  
  # Store the coefficients for opposition supporters
  plotdata$x[i * 2] <- coef(lm_model_opposition)["treated"]
  
  # Store the confidence intervals for incumbent supporters
  conf_int_90_incumbent <- confint(lm_model_incumbent, level = confidence_levels[1])
  plotdata$xmin90[i * 2 - 1] <- conf_int_90_incumbent["treated", 1]
  plotdata$xmax90[i * 2 - 1] <- conf_int_90_incumbent["treated", 2]
  
  conf_int_95_incumbent <- confint(lm_model_incumbent, level = confidence_levels[2])
  plotdata$xmin95[i * 2 - 1] <- conf_int_95_incumbent["treated", 1]
  plotdata$xmax95[i * 2 - 1] <- conf_int_95_incumbent["treated", 2]
  
  # Store the confidence intervals for opposition supporters
  conf_int_90_opposition <- confint(lm_model_opposition, level = confidence_levels[1])
  plotdata$xmin90[i * 2] <- conf_int_90_opposition["treated", 1]
  plotdata$xmax90[i * 2] <- conf_int_90_opposition["treated", 2]
  
  conf_int_95_opposition <- confint(lm_model_opposition, level = confidence_levels[2])
  plotdata$xmin95[i * 2] <- conf_int_95_opposition["treated", 1]
  plotdata$xmax95[i * 2] <- conf_int_95_opposition["treated", 2]
}

coeffplot <- ggplot(plotdata,aes(x = label, y = x, color = Group, shape = Group)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_linerange(aes(ymin = xmin90,
                     ymax = xmax90),
                 linewidth = 1,
                 position = position_dodge(0.4)) +
  geom_linerange(aes(ymin = xmin95,
                     ymax = xmax95),
                 linewidth = 0.45,
                 position = position_dodge(0.4)) + 
  geom_point(size=3,position = position_dodge(0.4)) + 
  xlab("") +
  ylab("Effect of Announcement (LATE)") +
  theme(legend.title=element_blank()) +
  coord_flip() + 
  theme(axis.text.x = element_text(size = 8, angle = 45,hjust = 1, vjust = 1)) + 
  scale_color_manual(values = plot_colors)  +
  theme(legend.position = "top",
        text = element_text(size = 15),
        axis.text = element_text(size = 14))

pdf(file = "fig_c4a.pdf",width = 7.5,height = 6,onefile=F)

coeffplot

dev.off()

# Figure C.4b -------------------------------------------------------------

# Function to run rdrobust and retrieve results
run_rd_robust <- function(data, data_incumbent, data_opposition, running, running_incumbent, running_opposition) {
  rdrb_full_95 <- rdrobust(data, running, c = 0, masspoints = "adjust", all = TRUE, level = 95)
  rdrb_full_90 <- rdrobust(data, running, c = 0, masspoints = "adjust", all = TRUE, level = 90)
  rdrb_in_95 <- rdrobust(data_incumbent, running_incumbent, c = 0, masspoints = "adjust", all = TRUE, level = 95)
  rdrb_in_90 <- rdrobust(data_incumbent, running_incumbent, c = 0, masspoints = "adjust", all = TRUE, level = 90)
  rdrb_op_95 <- rdrobust(data_opposition, running_opposition, c = 0, masspoints = "adjust", all = TRUE, level = 95)
  rdrb_op_90 <- rdrobust(data_opposition, running_opposition, c = 0, masspoints = "adjust", all = TRUE, level = 90)
  result <- list(rdrb_full_95, rdrb_full_90, rdrb_in_95, rdrb_in_90, rdrb_op_95, rdrb_op_90)
  return(result)
}

# List of dependent variables
dependent_vars <- c("Q49","Q59","Q142","Q147","Q170","Q199","Q222")

# Run rdrobust for each dependent variable
rd_results <- lapply(
  dependent_vars,
  function(dep_var) {
    run_rd_robust(
      data[[dep_var]], data_incumbent[[dep_var]], data_opposition[[dep_var]],
      data$D_INTERVIEW, data_incumbent$D_INTERVIEW, data_opposition$D_INTERVIEW
    )
  }
)

# Initialize the dataframe
plotdata_robust <- data.frame(
  y = rep(1:7, each = 2),
  Group = rep(c("Incumbent supporters", "Opposition supporters")),
  x = rep(1:7, each = 2),
  xmin95 = NA,
  xmax95 = NA,
  xmin90 = NA,
  xmax90 = NA,
  question = rep(c(dependent_vars), each = 2),
  label = rep(c(
    "Satisfaction with life",
    "Trust in neighborhood",
    "Fear of losing job",
    "Worried about terrorism",
    "Only own religion acceptable",
    "Interest in politics",
    "Frequency of voting"
  ), each = 2)
)

# Populate the dataframe with coefficients and confidence intervals
for (i in seq_along(dependent_vars)) {
  plotdata_robust$x[i * 2 - 1] <- rd_results[[i]][[3]]$coef[3]
  plotdata_robust$x[i * 2] <- rd_results[[i]][[5]]$coef[3]
  
  plotdata_robust$xmin90[i * 2 - 1] <- rd_results[[i]][[4]]$ci["Robust", "CI Lower"]
  plotdata_robust$xmax90[i * 2 - 1] <- rd_results[[i]][[4]]$ci["Robust", "CI Upper"]
  plotdata_robust$xmin95[i * 2 - 1] <- rd_results[[i]][[3]]$ci["Robust", "CI Lower"]
  plotdata_robust$xmax95[i * 2 - 1] <- rd_results[[i]][[3]]$ci["Robust", "CI Upper"]
  
  plotdata_robust$xmin90[i * 2] <- rd_results[[i]][[6]]$ci["Robust", "CI Lower"]
  plotdata_robust$xmax90[i * 2] <- rd_results[[i]][[6]]$ci["Robust", "CI Upper"]
  plotdata_robust$xmin95[i * 2] <- rd_results[[i]][[5]]$ci["Robust", "CI Lower"]
  plotdata_robust$xmax95[i * 2] <- rd_results[[i]][[5]]$ci["Robust", "CI Upper"]
}

coeffplot <- ggplot(plotdata_robust,aes(x = label, y = x, color = Group, shape = Group)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_linerange(aes(ymin = xmin90,
                     ymax = xmax90),
                 linewidth = 1,
                 position = position_dodge(0.4)) +
  geom_linerange(aes(ymin = xmin95,
                     ymax = xmax95),
                 linewidth = 0.45,
                 position = position_dodge(0.4)) + 
  geom_point(size=3,position = position_dodge(0.4)) + 
  xlab("") +
  ylab("Effect of Announcement (LATE)") +
  theme(legend.title=element_blank()) +
  coord_flip() + 
  theme(axis.text.x = element_text(size = 8, angle = 45,hjust = 1, vjust = 1)) + 
  scale_color_manual(values = plot_colors) +
  theme(legend.position = "top",
        text = element_text(size = 15),
        axis.text = element_text(size = 14))

pdf(file = "fig_c4b.pdf",width = 7.5,height = 6,onefile=F)

coeffplot

dev.off()

# Table C.16 --------------------------------------------------------------

BWdata <- data_incumbent[data_incumbent$D_INTERVIEW<=4&data_incumbent$D_INTERVIEW>=-4,]

tt1 <- t.test(BWdata$Q262 ~ BWdata$treated) # age
tt2 <- t.test(BWdata$Q260==1 ~ BWdata$treated) # gender
tt3 <- t.test(BWdata$H_URBRURAL==1 ~ BWdata$treated) # urban/rural
tt4 <- t.test(BWdata$Q279 ~ BWdata$treated) # employment status
tt5 <- t.test(BWdata$Q202 ~ BWdata$treated) # TV consumption

table <- map_df(list(tt1, tt2, tt3, tt4, tt5), tidy)
balance_incumbent <- table[c("estimate", "statistic", "p.value", "conf.low", "conf.high")]
balance_incumbent$variable <- c("Age","Gender","Urban/rural","Employment status","TV consumption")
write.csv(balance_incumbent, file = "tab_c16.csv",row.names = F)

# Table C.17 --------------------------------------------------------------

BWdata <- data_opposition[data_opposition$D_INTERVIEW<=4&data_opposition$D_INTERVIEW>=-4,]

tt1 <- t.test(BWdata$Q262 ~ BWdata$treated) # age
tt2 <- t.test(BWdata$Q260==1 ~ BWdata$treated) # gender
tt3 <- t.test(BWdata$H_URBRURAL==1 ~ BWdata$treated) # urban/rural
tt4 <- t.test(BWdata$Q279 ~ BWdata$treated) # employment status
tt5 <- t.test(BWdata$Q202 ~ BWdata$treated) # TV consumption

table <- map_df(list(tt1, tt2, tt3, tt4, tt5), tidy)
balance_opposition <- table[c("estimate", "statistic", "p.value", "conf.low", "conf.high")]
balance_opposition$variable <- c("Age","Gender","Urban/rural","Employment status","TV consumption")
write.csv(balance_opposition, file = "tab_c17.csv",row.names = F)

# Tables C.18 - C.28 ------------------------------------------------------

create_lm <- function(data, data_incumbent, data_opposition, dep_var) {
  lm1 <- lm(as.formula(paste0(dep_var, " ~ treated * D_INTERVIEW + Q262 + Q260 + H_URBRURAL + Q279 + Q202")), data = data[data$D_INTERVIEW >= -4 & data$D_INTERVIEW <= 4, ])
  lm2 <- lm(as.formula(paste0(dep_var, " ~ treated * D_INTERVIEW + Q262 + Q260 + H_URBRURAL + Q279 + Q202")), data = data_incumbent[data_incumbent$D_INTERVIEW >= -4 & data_incumbent$D_INTERVIEW <= 4, ])
  lm3 <- lm(as.formula(paste0(dep_var, " ~ treated * D_INTERVIEW + Q262 + Q260 + H_URBRURAL + Q279 + Q202")), data = data_opposition[data_opposition$D_INTERVIEW >= -5 & data_opposition$D_INTERVIEW <= 5, ])
  return(list(lm1, lm2, lm3))
}

# List of dependent variables
dependent_vars <- paste0("Q", 224:234)

# List to store the linear models
lm_list <- lapply(dependent_vars, function(dep_var) create_lm(data, data_incumbent, data_opposition, dep_var))

# Define the confidence levels
confidence_levels <- c(0.90, 0.95)

# Initialize the dataframe
plotdata <- data.frame(
  y = rep(1:11, each = 2),
  Group = rep(c("Incumbent supporters", "Opposition supporters"), times = 11),
  x = rep(1:11, each = 2),
  xmin95 = NA,
  xmax95 = NA,
  xmin90 = NA,
  xmax90 = NA,
  question = rep(paste0("Q", 224:234), each = 2),
  label = rep(c(
    "Votes counted fairly",
    "Opposition prevented from running",
    "TV favors governing party",
    "Voters bribed",
    "Journalists provide fair coverage",
    "Election officials fair",
    "Rich people buy elections",
    "Voters threatened with violence",
    "Voters offered genuine choice",
    "Women have equal opportunities to run",
    "Honest elections make a difference"
  ), each = 2)
)

for (i in seq_along(dependent_vars)) {
  dep_var <- dependent_vars[i]
  lm_model_incumbent <- lm_list[[i]][[2]]  # Access the second linear model
  lm_model_opposition <- lm_list[[i]][[3]]  # Access the third linear model
  
  # Store the coefficients for incumbent supporters
  plotdata$x[i * 2 - 1] <- coef(lm_model_incumbent)["treated"]
  
  # Store the coefficients for opposition supporters
  plotdata$x[i * 2] <- coef(lm_model_opposition)["treated"]
  
  # Store the confidence intervals for incumbent supporters
  conf_int_90_incumbent <- confint(lm_model_incumbent, level = confidence_levels[1])
  plotdata$xmin90[i * 2 - 1] <- conf_int_90_incumbent["treated", 1]
  plotdata$xmax90[i * 2 - 1] <- conf_int_90_incumbent["treated", 2]
  
  conf_int_95_incumbent <- confint(lm_model_incumbent, level = confidence_levels[2])
  plotdata$xmin95[i * 2 - 1] <- conf_int_95_incumbent["treated", 1]
  plotdata$xmax95[i * 2 - 1] <- conf_int_95_incumbent["treated", 2]
  
  # Store the confidence intervals for opposition supporters
  conf_int_90_opposition <- confint(lm_model_opposition, level = confidence_levels[1])
  plotdata$xmin90[i * 2] <- conf_int_90_opposition["treated", 1]
  plotdata$xmax90[i * 2] <- conf_int_90_opposition["treated", 2]
  
  conf_int_95_opposition <- confint(lm_model_opposition, level = confidence_levels[2])
  plotdata$xmin95[i * 2] <- conf_int_95_opposition["treated", 1]
  plotdata$xmax95[i * 2] <- conf_int_95_opposition["treated", 2]
}

# Define output filenames
output_files <- character(11)
for (i in 18:28) {
  output_files[i-17] <- paste0("tab_c", i, ".tex")
}

# Loop over lm_list and generate corresponding tables
for (i in seq_along(lm_list)) {
  texreg(lm_list[[i]],
         file = output_files[i], # Specify output filename
         include.ci = FALSE,
         digits = 3,
         custom.coef.map = list(
           "(Intercept)" = "Intercept",
           "treated" = "Treated",
           "D_INTERVIEW" = "Running variable",
           "treated:D_INTERVIEW" = "Treated $\\times$ running variable",
           "Q262" = "Age",
           "Q260" = "Gender",
           "H_URBRURAL" = "Urban/rural",
           "Q279" = "Employment status",
           "Q202" = "TV consumption"
         ),
         custom.header = list("Pooled" = 1,
                              "Incumbent" = 2,
                              "Opposition" = 3),
         booktabs = TRUE,
         use.packages = FALSE,
         include.rmse = FALSE
  )
}

# Table C.29 --------------------------------------------------------------

BWdata <- data_incumbent[data_incumbent$D_INTERVIEW<=4&data_incumbent$D_INTERVIEW>=-4,]

tt1 <- t.test(as.numeric(!is.na(BWdata$Q224))~ BWdata$treated)
tt2 <- t.test(as.numeric(!is.na(BWdata$Q225))~ BWdata$treated)
tt3 <- t.test(as.numeric(!is.na(BWdata$Q226))~ BWdata$treated)
tt4 <- t.test(as.numeric(!is.na(BWdata$Q227))~ BWdata$treated)
tt5 <- t.test(as.numeric(!is.na(BWdata$Q228))~ BWdata$treated)
tt6 <- t.test(as.numeric(!is.na(BWdata$Q229))~ BWdata$treated)
tt7 <- t.test(as.numeric(!is.na(BWdata$Q230))~ BWdata$treated)
tt8 <- t.test(as.numeric(!is.na(BWdata$Q231))~ BWdata$treated)
tt9 <- t.test(as.numeric(!is.na(BWdata$Q232))~ BWdata$treated)
tt10 <- t.test(as.numeric(!is.na(BWdata$Q233))~ BWdata$treated)
tt11 <- t.test(as.numeric(!is.na(BWdata$Q234))~ BWdata$treated)

table <- map_df(list(tt1, tt2, tt3, tt4, tt5, tt6, tt7, tt8, tt9, tt10, tt11), tidy)
balance_na_inucumbent <- table[c("estimate", "statistic", "p.value", "conf.low", "conf.high")]

balance_na_inucumbent$variable <- c(
  "Votes counted fairly",
  "Opposition prevented from running",
  "TV favors governing party",
  "Voters bribed",
  "Journalists provide fair coverage",
  "Election officials fair",
  "Rich people buy elections",
  "Voters threatened with violence",
  "Voters offered genuine choice",
  "Women have equal opportunities to run",
  "Honest elections make a difference"
)

write.csv(balance_na_inucumbent, file = "tab_c29.csv")

# Table C.30 --------------------------------------------------------------

BWdata <- data_opposition[data_opposition$D_INTERVIEW<=4&data_opposition$D_INTERVIEW>=-4,]
tt1 <- t.test(as.numeric(!is.na(BWdata$Q224))~ BWdata$treated)
tt2 <- t.test(as.numeric(!is.na(BWdata$Q225))~ BWdata$treated)
tt3 <- t.test(as.numeric(!is.na(BWdata$Q226))~ BWdata$treated)
tt4 <- t.test(as.numeric(!is.na(BWdata$Q227))~ BWdata$treated)
tt5 <- t.test(as.numeric(!is.na(BWdata$Q228))~ BWdata$treated)
tt6 <- t.test(as.numeric(!is.na(BWdata$Q229))~ BWdata$treated)
tt7 <- t.test(as.numeric(!is.na(BWdata$Q230))~ BWdata$treated)
tt8 <- t.test(as.numeric(!is.na(BWdata$Q231))~ BWdata$treated)
tt9 <- t.test(as.numeric(!is.na(BWdata$Q232))~ BWdata$treated)
tt10 <- t.test(as.numeric(!is.na(BWdata$Q233))~ BWdata$treated)
tt11 <- t.test(as.numeric(!is.na(BWdata$Q234))~ BWdata$treated)

table <- map_df(list(tt1, tt2, tt3, tt4, tt5, tt6, tt7, tt8, tt9, tt10, tt11), tidy)
balance_na_opposition <- table[c("estimate", "statistic", "p.value", "conf.low", "conf.high")]
balance_na_opposition$variable <- c(
  "Votes counted fairly",
  "Opposition prevented from running",
  "TV favors governing party",
  "Voters bribed",
  "Journalists provide fair coverage",
  "Election officials fair",
  "Rich people buy elections",
  "Voters threatened with violence",
  "Voters offered genuine choice",
  "Women have equal opportunities to run",
  "Honest elections make a difference"
)

write.csv(balance_na_opposition, file = "tab_c30.csv")

# Figure C.5 --------------------------------------------------------------

Y = data_incumbent$Q232
R = data_incumbent$running

pdf(file = "fig_c5.pdf",
    width = 4,
    height = 4,
    onefile=F)

tmp <- rdpower(data=cbind(Y,R), masspoints = "adjust",alpha=0.1,plot=T,
               samph=c(10,10))

dev.off()

# Figure C.6 --------------------------------------------------------------

Y = data_opposition$Q232
R = data_opposition$running

pdf(file = "fig_c6.pdf",
    width = 4,
    height = 4,
    onefile=F)

tmp <- rdpower(data=cbind(Y,R), masspoints = "adjust",alpha=0.1,plot=T,
               samph=c(10,10))

dev.off()

# Figure C.7a -------------------------------------------------------------

create_lm <- function(data, data_incumbent, data_opposition, dep_var) {
  lm1 <- lm(as.formula(paste0(dep_var, " ~ treated * D_INTERVIEW")), data = data[data$D_INTERVIEW >= -10 & data$D_INTERVIEW <= 10, ])
  lm2 <- lm(as.formula(paste0(dep_var, " ~ treated * D_INTERVIEW")), data = data_incumbent[data_incumbent$D_INTERVIEW >= -10 & data_incumbent$D_INTERVIEW <= 10, ])
  lm3 <- lm(as.formula(paste0(dep_var, " ~ treated * D_INTERVIEW")), data = data_opposition[data_opposition$D_INTERVIEW >= -10 & data_opposition$D_INTERVIEW <= 10, ])
  return(list(lm1, lm2, lm3))
}

# List of dependent variables
dependent_vars <- paste0("Q", 224:234)

# List to store the linear models
lm_list <- lapply(dependent_vars, function(dep_var) create_lm(data, data_incumbent, data_opposition, dep_var))

# Define the confidence levels
confidence_levels <- c(0.90, 0.95)

# Initialize the dataframe
plotdata <- data.frame(
  y = rep(1:11, each = 2),
  Group = rep(c("Incumbent supporters", "Opposition supporters"), times = 11),
  x = rep(1:11, each = 2),
  xmin95 = NA,
  xmax95 = NA,
  xmin90 = NA,
  xmax90 = NA,
  question = rep(paste0("Q", 224:234), each = 2),
  label = rep(c(
    "Votes counted fairly",
    "Opposition prevented from running",
    "TV favors governing party",
    "Voters bribed",
    "Journalists provide fair coverage",
    "Election officials fair",
    "Rich people buy elections",
    "Voters threatened with violence",
    "Voters offered genuine choice",
    "Women have equal opportunities to run",
    "Honest elections make a difference"
  ), each = 2)
)

for (i in seq_along(dependent_vars)) {
  dep_var <- dependent_vars[i]
  lm_model_incumbent <- lm_list[[i]][[2]]  # Access the second linear model
  lm_model_opposition <- lm_list[[i]][[3]]  # Access the third linear model
  
  # Store the coefficients for incumbent supporters
  plotdata$x[i * 2 - 1] <- coef(lm_model_incumbent)["treated"]
  
  # Store the coefficients for opposition supporters
  plotdata$x[i * 2] <- coef(lm_model_opposition)["treated"]
  
  # Store the confidence intervals for incumbent supporters
  conf_int_90_incumbent <- confint(lm_model_incumbent, level = confidence_levels[1])
  plotdata$xmin90[i * 2 - 1] <- conf_int_90_incumbent["treated", 1]
  plotdata$xmax90[i * 2 - 1] <- conf_int_90_incumbent["treated", 2]
  
  conf_int_95_incumbent <- confint(lm_model_incumbent, level = confidence_levels[2])
  plotdata$xmin95[i * 2 - 1] <- conf_int_95_incumbent["treated", 1]
  plotdata$xmax95[i * 2 - 1] <- conf_int_95_incumbent["treated", 2]
  
  # Store the confidence intervals for opposition supporters
  conf_int_90_opposition <- confint(lm_model_opposition, level = confidence_levels[1])
  plotdata$xmin90[i * 2] <- conf_int_90_opposition["treated", 1]
  plotdata$xmax90[i * 2] <- conf_int_90_opposition["treated", 2]
  
  conf_int_95_opposition <- confint(lm_model_opposition, level = confidence_levels[2])
  plotdata$xmin95[i * 2] <- conf_int_95_opposition["treated", 1]
  plotdata$xmax95[i * 2] <- conf_int_95_opposition["treated", 2]
}

coeffplot <- ggplot(plotdata,aes(x = label, y = x, color = Group, shape = Group)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_linerange(aes(ymin = xmin90,
                     ymax = xmax90),
                 linewidth = 1,
                 position = position_dodge(0.4)) +
  geom_linerange(aes(ymin = xmin95,
                     ymax = xmax95),
                 linewidth = 0.45,
                 position = position_dodge(0.4)) + 
  geom_point(size=3,position = position_dodge(0.4)) + 
  xlab("") +
  ylab("Effect of Announcement (LATE)") +
  theme(legend.title=element_blank()) +
  coord_flip() + 
  theme(axis.text.x = element_text(size = 8, angle = 45,hjust = 1, vjust = 1)) + 
  scale_color_manual(values = plot_colors)  +
  theme(legend.position = "top",
        text = element_text(size = 15),
        axis.text = element_text(size = 14))

pdf(file = "fig_c7a.pdf",width = 7.5,height = 6,onefile=F)

coeffplot

dev.off()

# Figure C.7b -------------------------------------------------------------

# Function to run rdrobust and retrieve results
run_rd_robust <- function(data, data_incumbent, data_opposition, running, running_incumbent, running_opposition) {
  rdrb_full_95 <- rdrobust(data, running, c = 0, masspoints = "adjust", all = TRUE, level = 95, h = 10, b = 10)
  rdrb_full_90 <- rdrobust(data, running, c = 0, masspoints = "adjust", all = TRUE, level = 90, h = 10, b = 10)
  rdrb_in_95 <- rdrobust(data_incumbent, running_incumbent, c = 0, masspoints = "adjust", all = TRUE, level = 95, h = 10, b = 10)
  rdrb_in_90 <- rdrobust(data_incumbent, running_incumbent, c = 0, masspoints = "adjust", all = TRUE, level = 90, h = 10, b = 10)
  rdrb_op_95 <- rdrobust(data_opposition, running_opposition, c = 0, masspoints = "adjust", all = TRUE, level = 95, h = 10, b = 10)
  rdrb_op_90 <- rdrobust(data_opposition, running_opposition, c = 0, masspoints = "adjust", all = TRUE, level = 90, h = 10, b = 10)
  result <- list(rdrb_full_95, rdrb_full_90, rdrb_in_95, rdrb_in_90, rdrb_op_95, rdrb_op_90)
  return(result)
}

# List of dependent variables
dependent_vars <- paste0("Q", 224:234)

# Run rdrobust for each dependent variable
rd_results <- lapply(
  dependent_vars,
  function(dep_var) {
    run_rd_robust(
      data[[dep_var]], data_incumbent[[dep_var]], data_opposition[[dep_var]],
      data$D_INTERVIEW, data_incumbent$D_INTERVIEW, data_opposition$D_INTERVIEW
    )
  }
)

# Initialize the dataframe
plotdata_robust <- data.frame(
  y = rep(1:11, each = 2),
  Group = rep(c("Incumbent supporters", "Opposition supporters"), times = 11),
  x = rep(1:11, each = 2),
  xmin95 = NA,
  xmax95 = NA,
  xmin90 = NA,
  xmax90 = NA,
  question = rep(paste0("Q", 224:234), each = 2),
  label = rep(c(
    "Votes counted fairly",
    "Opposition prevented from running",
    "TV favors governing party",
    "Voters bribed",
    "Journalists provide fair coverage",
    "Election officials fair",
    "Rich people buy elections",
    "Voters threatened with violence",
    "Voters offered genuine choice",
    "Women have equal opportunities to run",
    "Honest elections make a difference"
  ), each = 2)
)

# Populate the dataframe with coefficients and confidence intervals
for (i in seq_along(dependent_vars)) {
  plotdata_robust$x[i * 2 - 1] <- rd_results[[i]][[3]]$coef[3]
  plotdata_robust$x[i * 2] <- rd_results[[i]][[5]]$coef[3]
  
  plotdata_robust$xmin90[i * 2 - 1] <- rd_results[[i]][[4]]$ci["Robust", "CI Lower"]
  plotdata_robust$xmax90[i * 2 - 1] <- rd_results[[i]][[4]]$ci["Robust", "CI Upper"]
  plotdata_robust$xmin95[i * 2 - 1] <- rd_results[[i]][[3]]$ci["Robust", "CI Lower"]
  plotdata_robust$xmax95[i * 2 - 1] <- rd_results[[i]][[3]]$ci["Robust", "CI Upper"]
  
  plotdata_robust$xmin90[i * 2] <- rd_results[[i]][[6]]$ci["Robust", "CI Lower"]
  plotdata_robust$xmax90[i * 2] <- rd_results[[i]][[6]]$ci["Robust", "CI Upper"]
  plotdata_robust$xmin95[i * 2] <- rd_results[[i]][[5]]$ci["Robust", "CI Lower"]
  plotdata_robust$xmax95[i * 2] <- rd_results[[i]][[5]]$ci["Robust", "CI Upper"]
}

coeffplot <- ggplot(plotdata_robust,aes(x = label, y = x, color = Group, shape = Group)) +
  geom_hline(yintercept = 0, 
             colour = gray(1/2), lty = 2) +
  geom_linerange(aes(ymin = xmin90,
                     ymax = xmax90),
                 linewidth = 1,
                 position = position_dodge(0.4)) +
  geom_linerange(aes(ymin = xmin95,
                     ymax = xmax95),
                 linewidth = 0.45,
                 position = position_dodge(0.4)) + 
  geom_point(size=3,position = position_dodge(0.4)) + 
  xlab("") +
  ylab("Effect of Announcement (LATE)") +
  theme(legend.title=element_blank()) +
  coord_flip() + 
  theme(axis.text.x = element_text(size = 8, angle = 45,hjust = 1, vjust = 1)) + 
  scale_color_manual(values = plot_colors) +
  theme(legend.position = "top",
        text = element_text(size = 15),
        axis.text = element_text(size = 14))

coeffplot

pdf(file = "fig_c7b.pdf",width = 7.5,height = 6,onefile=F)

coeffplot

dev.off()